*This code calculates horizontal shifts and analyzes whether it is an admissable shift (i.e. both monotonic relative Engel curves)
*this file allows for a bootstrap (bs input) but just takes curves from engel_lpoly_loop_v1_bs.do



qui {
**********need to split up and include the loop

local start = `1'
local finish = `2'
local Gi_spec "`3'"
local share_spec "`4'"
local band_spec "`5'"
global dump "`6'"
global dropbox "`7'"
local iteration = `8'
local foldno = "`9'"
local bs = `10'





set more off
set matsize 11000


foreach DMI in   `Gi_spec' {  


cd "${dump}/Engel`foldno'/"    


********************************************************************************
*Step 2: 
*Pull in one market-round at a time from the dataset above. 
*There will be the number of observations=number of households observed in market (possibly less than 100).
********************************************************************************

foreach band in  wbw`band_spec'  {  
noi di "*****************************************************************************************"
noi di "`band'"
noi di "*****************************************************************************************"
forval market=`start'/`finish' {
noi di "`market' "





foreach share in `share_spec'   {  

capture confirm file Engel_rindexes_`band_spec'_`share'_`DMI'_`market'_bs`bs'.dta
if _rc!=0 {

cap use "Engel_`band_spec'_`share_spec'_`DMI'_`market'_bs`bs'.dta", clear
if _rc==0 {



egen avg_count_hh=rowmean(r??count_hh)


*clean up
renvars r??count_hh *`band'_`share'*, prefix(x)
cap drop r??count* 
cap drop r??*wbw*
renvars xr??count_hh x*`band'_`share'*, presub("x")



*first get rid of noisy estimates and interpolate if they are interior points
foreach rnd in 43  55 { 
foreach var of varlist   r`rnd'lp_w* {

local varse=regexr("`var'","lp_","se_")


***now if 5-95 are montonic, but above is not, fix it
local namer3=regexr("`var'","_`share'","d_`share'")
local namerse3=regexr("`namer3'","lp_","se_")
gen slope=(`var'[_n+1]>`var'[_n] & `var'[_n]>`var'[_n-1])-(`var'[_n+1]<`var'[_n] & `var'[_n]<`var'[_n-1]) if _n<=101 & `var'[_n]!=. & `var'[_n-1]!=. & `var'[_n+1]!=.
replace slope=(`var'[_n+1]>`var'[_n])-(`var'[_n+1]<`var'[_n]) if _n<=101 & `var'[_n-1]==. & `var'[_n]!=. & `var'[_n+1]!=.
replace slope=(`var'[_n]>`var'[_n-1])-( `var'[_n]<`var'[_n-1]) if _n<=101 & `var'[_n+1]==. & `var'[_n]!=. & `var'[_n-1]!=.

sum slope if _n<=101, mean
local monotonic_0199=(r(max)==1 & r(min)==1)-(r(max)==-1 & r(min)==-1)
sum slope if _n>5 & _n<97, mean
local monotonic_0595=(r(max)==1 & r(min)==1)-(r(max)==-1 & r(min)==-1)
*flag for full monotonicity (-1 negative, 1 positive)

if `monotonic_0595'!=0 & `monotonic_0199'==0 {
gen `namer3'=`var'  if slope[_n]==`monotonic_0595' | (slope[_n]==0 & slope[_n-1]==`monotonic_0595') | (slope[_n]==0 & slope[_n+1]==`monotonic_0595')
replace `namer3'=. if (`namer3'[_n+1]==. | `namer3'[_n+2]==. | `namer3'[_n+3]==. | `namer3'[_n+4]==.) & _n<=5
replace `namer3'=. if (`namer3'[_n-1]==. | `namer3'[_n-2]==. | `namer3'[_n-3]==. | `namer3'[_n-4]==.) & _n>=97 & _n<=101
*** this just drops the non monotonic bit at ends
}
else {
gen `namer3'=`var'
}
cap gen `namerse3'=`varse' if `namer3'!=.
drop slope


*always extrapolate at end but also clean up non monotonic bits before then
local namer6=regexr("`var'","_`share'","ce_`share'")

gen xtemp=`namer3' if  _n>3 & _n<99   
mipolate xtemp r`rnd'predict_lnmpcew, spline epolate gen(`namer6')
drop xtemp




}
}





*clean up
drop *d_`share'*
drop *`band_spec'_`share'*




*calculate nominal income change at each percentile
foreach pct in 10 20 30 40 50 60 70 80 90  {
foreach rnd in 43 55 {  
local income`rnd'`pct'=r`rnd'predict_lnmpcew[`pct'+1]
}

foreach rnd1 in 43 55 {  
foreach rnd2 in 43  55 {  
if `rnd1'!=`rnd2' {
gen W`rnd1'`rnd2'_lnmpcew_p`pct'=(`income`rnd2'`pct''-`income`rnd1'`pct'')
}
}
}

}





local newobs=_N+9
set obs `newobs'
replace percentile=999505 if _n==`newobs'-2
replace percentile=999703 if _n==`newobs'-1
replace percentile=999901 if _n==`newobs'
*whether curves are monotonic

replace percentile=99505 if _n==`newobs'-5
replace percentile=99703 if _n==`newobs'-4
replace percentile=99901 if _n==`newobs'-3
*reason if no cross

replace percentile=9505 if _n==`newobs'-8
replace percentile=9703 if _n==`newobs'-7
replace percentile=9901 if _n==`newobs'-6
*value if unique crossing







foreach var of varlist   r43lp_wbw*  { //
local var00=regexr("`var'","^r43","")

foreach rnd in 43  55 {  
local var`rnd'=regexr("`var'","^r43","r`rnd'")
local se`rnd'=regexr("`var`rnd''","lp","se")
*this is key code as used to judge if monotonic
gen slope`rnd'=(`var`rnd''[_n+1]>`var`rnd''[_n] & `var`rnd''[_n]>`var`rnd''[_n-1])-(`var`rnd''[_n+1]<`var`rnd''[_n] & `var`rnd''[_n]<`var`rnd''[_n-1]) if _n<=101 & `var`rnd''[_n]!=. & `var`rnd''[_n-1]!=. & `var`rnd''[_n+1]!=.
replace slope`rnd'=(`var`rnd''[_n+1]>`var`rnd''[_n])-(`var`rnd''[_n+1]<`var`rnd''[_n]) if _n<=101 & `var`rnd''[_n-1]==. & `var`rnd''[_n]!=. & `var`rnd''[_n+1]!=.
replace slope`rnd'=(`var`rnd''[_n]>`var`rnd''[_n-1])-( `var`rnd''[_n]<`var`rnd''[_n-1]) if _n<=101 & `var`rnd''[_n+1]==. & `var`rnd''[_n]!=. & `var`rnd''[_n-1]!=.
}


foreach rnd in 43  55 { 
sum slope`rnd' if _n<=101, mean
local monotonic_0199_`rnd'=(r(max)==1 & r(min)==1)+2*(r(max)==-1 & r(min)==-1)
*flag for full monotonicity (2 negative, 1 positive)
}


foreach pct in 10 20 30 40 50 60 70 80 90  {

foreach rnd in 43  55 { 
local bshare`rnd'`pct'=`var`rnd''[`pct'+1]
local income`rnd'`pct'=r`rnd'predict_lnmpcew[`pct'+1]
local slope`rnd'`pct'=slope`rnd'[`pct'+1]
cap local se`rnd'`pct'=`se`rnd''[`pct'+1]
}








****first go forward (e.g. 43-->55)
foreach rnd in  55 { 

gen P43`rnd'`var00'_p`pct'=(r`rnd'predict_lnmpcew[_n-1] + (r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1])*(`bshare43`pct''-`var`rnd''[_n-1])/(`var`rnd''-`var`rnd''[_n-1]) - `income43`pct'') ///
if  (`var`rnd''>`bshare43`pct'' & `var`rnd''[_n-1]<`bshare43`pct'' & `var`rnd''!=. & `var`rnd''[_n-1]!=. & `slope43`pct''==slope`rnd' & `slope43`pct''==1)  

replace P43`rnd'`var00'_p`pct'=(r`rnd'predict_lnmpcew[_n-1] + (r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1])*(`var`rnd''[_n-1]-`bshare43`pct'')/(`var`rnd''[_n-1]-`var`rnd'') - `income43`pct'') ///
if  (`var`rnd''[_n-1]>`bshare43`pct'' & `var`rnd''[_n]<`bshare43`pct'' & `var`rnd''!=. & `var`rnd''[_n-1]!=. & `slope43`pct''==slope`rnd' & `slope43`pct''==-1)   


*now get slopes at origin and dest
gen T43`rnd'`var00'_p`pct'=(`var`rnd''[_n]-`var`rnd''[_n-1])/(r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1]) if P43`rnd'`var00'_p`pct'!=.

gen byte S43`rnd'`var00'_p`pct'x=0.5*(`var43'[`pct'+2]-`var43'[`pct'+1])/(r43predict_lnmpcew[`pct'+2]-r43predict_lnmpcew[`pct'+1]) + 0.5*(`var43'[`pct'+1]-`var43'[`pct'])/(r43predict_lnmpcew[`pct'+1]-r43predict_lnmpcew[`pct']) 
replace S43`rnd'`var00'_p`pct'x=(`var43'[`pct'+2]-`var43'[`pct'+1])/(r43predict_lnmpcew[`pct'+2]-r43predict_lnmpcew[`pct'+1]) if S43`rnd'`var00'_p`pct'x==.
replace S43`rnd'`var00'_p`pct'x=(`var43'[`pct'+1]-`var43'[`pct'])/(r43predict_lnmpcew[`pct'+1]-r43predict_lnmpcew[`pct']) if S43`rnd'`var00'_p`pct'x==.



sum `var`rnd'' if _n<=101, mean
local xmax0199=r(max)
local xmin0199=r(min)



foreach var of varlist ?43`rnd'`var00'_p`pct' {

***********full range
sum `var'  if _n<=101, mean
if r(N)==1 {
replace `var'=r(max) if _n==_N-6
*this is unique crossing value
}
else if r(N)==0 {
	replace `var'= (`bshare43`pct''>`xmax0199' & `monotonic_0199_`rnd''==2) + (`bshare43`pct''<`xmin0199' & `monotonic_0199_`rnd''==1) ///
			-(`bshare43`pct''>`xmax0199' & `monotonic_0199_`rnd''==1) - (`bshare43`pct''<`xmin0199' & `monotonic_0199_`rnd''==2) if _n==_N-3
	replace `var'=99  if `var'==0 & _n==_N-3
	*reason for no cross (for 4355, value=1 if no real income match for poor, value=-1 if no real income match for rich, both require monotonicity in round 55)
}
else if r(N)>1 {
	replace `var'=0 if _n==_N-3
	*reason for no cross is multiple crossings
}
replace `var'=(43`monotonic_0199_43'`rnd'`monotonic_0199_`rnd'') if _n==_N
*monotonicity so 43X55X tells us about monotonicity, where X=0 mean fail, X=1 means positive slope, X=2 means negative slope



}


}





****now go backwards (e.g. 55-->43)




foreach rnd in 43 {
gen P55`rnd'`var00'_p`pct'=(r`rnd'predict_lnmpcew[_n-1] + (r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1])*(`bshare55`pct''-`var`rnd''[_n-1])/(`var`rnd''-`var`rnd''[_n-1]) - `income55`pct'') ///
if  (`var`rnd''>`bshare55`pct'' & `var`rnd''[_n-1]<`bshare55`pct'' & `var`rnd''!=. & `var`rnd''[_n-1]!=. & `slope55`pct''==slope`rnd' & `slope55`pct''==1)  

replace P55`rnd'`var00'_p`pct'=(r`rnd'predict_lnmpcew[_n-1] + (r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1])*(`var`rnd''[_n-1]-`bshare55`pct'')/(`var`rnd''[_n-1]-`var`rnd'') - `income55`pct'') ///
if  (`var`rnd''[_n-1]>`bshare55`pct'' & `var`rnd''[_n]<`bshare55`pct'' & `var`rnd''!=. & `var`rnd''[_n-1]!=. & `slope55`pct''==slope`rnd' & `slope55`pct''==-1)   

gen T55`rnd'`var00'_p`pct'=(`var`rnd''[_n]-`var`rnd''[_n-1])/(r`rnd'predict_lnmpcew[_n]-r`rnd'predict_lnmpcew[_n-1]) if P55`rnd'`var00'_p`pct'!=.

gen byte S55`rnd'`var00'_p`pct'x=0.5*(`var55'[`pct'+2]-`var55'[`pct'+1])/(r55predict_lnmpcew[`pct'+2]-r55predict_lnmpcew[`pct'+1]) + 0.5*(`var55'[`pct'+1]-`var55'[`pct'])/(r55predict_lnmpcew[`pct'+1]-r55predict_lnmpcew[`pct']) 
replace S55`rnd'`var00'_p`pct'x=(`var55'[`pct'+2]-`var55'[`pct'+1])/(r55predict_lnmpcew[`pct'+2]-r55predict_lnmpcew[`pct'+1]) if S55`rnd'`var00'_p`pct'x==.
replace S55`rnd'`var00'_p`pct'x=(`var55'[`pct'+1]-`var55'[`pct'])/(r55predict_lnmpcew[`pct'+1]-r55predict_lnmpcew[`pct']) if S55`rnd'`var00'_p`pct'x==.


sum `var`rnd'' if _n<=101, mean
local xmax0199=r(max)
local xmin0199=r(min)




foreach var of varlist ?55`rnd'`var00'_p`pct' {

***********full range
sum `var'  if _n<=101, mean
if r(N)==1 {
replace `var'=r(max) if _n==_N-6
*this is unique crossing value
}
else if r(N)==0 {
	replace `var'= (`bshare55`pct''>`xmax0199' & `monotonic_0199_`rnd''==2) + (`bshare55`pct''<`xmin0199' & `monotonic_0199_`rnd''==1) ///
			-(`bshare55`pct''>`xmax0199' & `monotonic_0199_`rnd''==1) - (`bshare55`pct''<`xmin0199' & `monotonic_0199_`rnd''==2) if _n==_N-3
	replace `var'=99  if `var'==0 & _n==_N-3
	*reason for no cross (for 5555, value=1 if no real income match for poor, value=-1 if no real income match for rich, both require monotonicity in round 55)
}
else if r(N)>1 {
	replace `var'=0 if _n==_N-3
	*reason for no cross is multiple crossings
}
replace `var'=(55`monotonic_0199_55'`rnd'`monotonic_0199_`rnd'') if _n==_N
*monotonicity so 55X55X tells us about monotonicity, where X=0 mean fail, X=1 means positive slope, X=2 means negative slope



}

}




}




foreach rnd in 43  55 {  
cap drop `var`rnd'' slope`rnd'
cap drop `se`rnd'' 
}


}



foreach var of varlist market_id sector state43 district43 avg_count_hh r??total_wt r??count_hh W* {
replace `var'=`var'[_n-1] if `var'==.
}


drop if percentile<=100
compress
save Engel_rindexes_`band_spec'_`share'_`DMI'_`market'_bs`bs'.dta, replace

 

}


}

}
}

}

clear
set obs 100
gen test=1
save InstanceB_`band_spec'_`share_spec'_`DMI'_`iteration'_bs`bs'.dta, replace

} // end of qui


}

exit, clear STATA




